TD Dynamique adaptative - MDet

Modele ecologique

On considere un modèle de competition. Dans notre modèle, 21 phénotypes dont le trait X est distribué sur un intervalle allant de 0 à 2, sont en compétition. Chaque population de chaque phénotype possède un taux de croissance intrasèque identique r. La densité maximale est spécifique à chaque phénotype et dépendant uniquement de la valeur de son trait selon la formule :

\[ K(x) = max(K_0 - λ(x-x_0)^2, 0) \] L’intensité de la compétition entre deux phénotypes i et j vaut : \[ a(x_i, x_j) =exp(-\frac{(x_i-xj)^2}{2σ^2}) \]

La dynamique de chaque population s’écrit alors :

\[ \frac{dn_i}{dt} = r.n_i.(\sum_{j=1}^{n} \frac{a(x_i, x_j).n_j}{K(x_i)} ) \] Pour l’ensemble des manipulations, on fixe certaines valeurs par défaut aux paramètres.

#Parametres
i = 21 # Nombre d'especes
Xmin = 0
Xmax = 2
X0 = (Xmax - Xmin)/2
K0 = 1
lambda = 1
sigma = 1 #etalement de la competition
r = 1

On implemente les fonctions méthémqtiques du modèle.

############ Fonctions #############

K_cap = function(x, K0_ = K0, lambda_= lambda, x0 = X0){
   (K0_ - lambda_*(x-x0)**2)
}

a = function(x1,x2, sigma_a = sigma){ 
  exp( (-1/2*(x2-x1)**2)/sigma_a**2 )
}

fitness = function(x1, x2, r_ = r ){
  r_*(1-a(x1, x2)*(K_cap(x1)/K_cap(x2)))
}

Fitness d’invasion

Afin de comprendre ce qu’implique notre modèle avec les valeurs de ses paramètres initiaux, nous representaons graphiquement la valeur de la fitness relative dans le cas de l’appration d’un nouveau phénotype. Ces figure en trois dimensions fournissent le paysage adaptatif de notre modèle.

## Set range and domain of plot
x1_interval  <- seq(0.1, 1.9, length.out = 25);
x2_interval  <- seq(0.1, 1.9, length.out = 25);

## Interpolate surface

z  <- outer(x1_interval,x2_interval,
            FUN = fitness)

p  <- persp(x1_interval,x2_interval,z, theta = 30, phi = 20,
            col = "lightblue", shade = 0.4, ticktype = "detailed")

plot_ly() %>% add_surface(x = ~x1_interval, y = ~x2_interval, z = ~z)

Le Pairwise Invasibily Plot (PIP) fournit un résumé de la situation en représentant seulement le signe de la différence des fitness.

## Set range and domain of plot
x1_interval  <- seq(0.1, 1.9, length.out = 500);
x2_interval  <- seq(0.1, 1.9, length.out = 500);

## Interpolate surface

z  <- outer(x1_interval,x2_interval,
            FUN = fitness)

couleurs <- ifelse(z > 0, "lightblue", "black")

matrice_couleurs = couleurs
df <- data.frame(x = rep(1:nrow(matrice_couleurs), each = ncol(matrice_couleurs)),
                 y = rep(1:ncol(matrice_couleurs), times = nrow(matrice_couleurs)),
                 couleur = as.vector(matrice_couleurs))

ggplot(df, aes(x = x, y = y, fill = couleur)) +
  geom_tile() +
  scale_fill_identity() +
  labs(title = "PIP") +
  theme_minimal()+
  labs(fill = "Légende\nNoir = Négatif\nBleu = Positif")

Simulations

Pour simuler ces compétitions, le modèle est implémenté sous forme matricielle. La dynamique est simulé pour un temps de 1000 et pour deux valeurs de sigma.

############### Simuler LV pour un nombre arbitraire d'espece ############

time_max = 1000

#Conditions initiales
x = seq(from = 0, to = 2, length.out = i) # Valeurs des phenotypes
k = c()
for (i in 1:length(x)){
  k = c(c(k),c(K_cap(x[i])))
}
 # K de chaque espece
n = rep(1/i,i) #densite de pop initiale

# On met sous forme de matrice
N0 = t(as.matrix(n))
K = matrix(rep(k,i), ncol = i, nrow = i, byrow = FALSE)
X = matrix(rep(x,i), ncol = i, nrow = i, byrow = TRUE)

#Calcul de la matrice M
D = t(X) - X


#Resolution du systeme d'ODE
model <- function(t,N,sigma){
  
  M = exp(-0.5*D^2/(sigma^2))/t(K)
  
  dN <- r * N * (1 - N %*% M)
  
  return(list(dN))
}

ode <-ode(N0,0:time_max,model,parms = sigma)

sigma = 5
# Representation
#barplot(ode[100,1:i+1])
ode <- ode(N0,0:time_max,model,parms = sigma)

ode <- as.data.frame(ode)
#colnames(ode) <- c("t",'x1','x2','x3')
ode_plot = ode%>%
  pivot_longer(-time, values_to = "dens", names_to = "sp_ID")

ggplot(ode_plot)+
  geom_line(aes(time, dens, col = sp_ID))

sigma = 0.3
# Representation
#barplot(ode[100,1:i+1])
ode <-ode(N0,0:time_max,model,parms = sigma)

ode <- as.data.frame(ode)
#colnames(ode) <- c("t",'x1','x2','x3')
ode_plot = ode%>%
  pivot_longer(-time, values_to = "dens", names_to = "sp_ID")

ggplot(ode_plot)+
  geom_line(aes(time, dens, col = sp_ID))

On observe que le nombre de phénotypes sui persistent dans le temps dépend notamment de la valeurs de sigma. Pour des faibles valeurs, un plus grand nombres de phénotypes se maintiennent.

ode_plot$sp_ID <- as.numeric(ode_plot$sp_ID)
ode_plot$trait <- x[as.numeric(ode_plot$sp_ID)]

ggplot(ode_plot)+
  geom_raster(aes(trait, time,  fill=dens**2))+
  scale_fill_gradient2(low = "white" ,
                       high = "red")

Le script et le graphique suivant permettent d’étudier le nombre de phénotypes se maintenant en fonction de la valeur de sigma. Les valeurs de sigma testées sont comprises entre 0.1 et 2.

seuil <- 0.001
sigma_bank = seq(0.1,2,by = 0.05)

nb_especes <- c()

for (sigma_i in sigma_bank){{
  
  ode2 <-ode(N0,0:time_max,model,parms = sigma_i)
  nb = sum(ode2[nrow((ode2)),-1] > seuil)
  }
  nb_especes <- c(c(nb_especes),c(nb))
  
}

res = tibble(sigma_bank,nb_especes)

ggplot()+
  geom_point(data = res, aes(x = sigma_bank, y = nb_especes))+
  labs(title = " ",
       x = "Sigma",
       y = "Nombre d'espèces en régime stationnaire")

Modelisation de mutations

Dans cette partie,

######## Implementation de mutation ##########

#parametres

p = 0.1 #taux de la pop qui mute
T = 10000
t = 200
n = rep(0,i)
n[i/6] <- 1 # densites de pop initiale (une seule espece)
N0 <- t(as.matrix(n))

sigma = 0.25


#Initialisation des valeures de K
k = c()
for (i in 1:length(x)){
  k = c(c(k),c(K_cap(x[i])))
}
K = matrix(rep(k,i), ncol = i, nrow = i, byrow = FALSE)


A = exp(-0.5*D^2/(sigma^2))
M = A/t(K)

ode <-ode(N0,1:t,model,parms = sigma)
ode_mutation <- ode

#trouver un moyen de faire muter tout le monde... (peut etre metre un seuil de densité min..)

for (z in 2:(T/t)){
  N0 <- t(as.matrix(tail(ode,1)[,1:i+1]))
  a = runif(1)
  index = which.max(N0)
  if (a<0.5){ #mutation vers la gauche
    N0[max(0,index-1)] <- p*N0[index]
    N0[index] <- (1-p)*N0[index]
  } else { #mutation vers la droite
    N0[min(i,index+1)] <- p*N0[index]
    N0[index] <- (1-p)*N0[index]
  }
  ode <-ode(N0,(z*t):(z*t+t),model,parms = sigma)
  
  ode_mutation <- rbind(ode_mutation,ode)
}

ode_mutation_plot = as.data.frame(ode_mutation)%>%
  pivot_longer(-time, values_to = "dens", names_to = "sp_ID")

ode_mutation_plot$sp_ID <- as.numeric(ode_mutation_plot$sp_ID)
ode_mutation_plot$trait <- x[as.numeric(ode_mutation_plot$sp_ID)]

ggplot(ode_mutation_plot)+
  geom_raster(aes(trait, time,  fill=dens**2))+
  scale_fill_gradient2(low = "white" ,
                       high = "red") +
  main_theme

## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.

# ggplot(ode_mutation_plot)+
#   geom_tile(aes(trait, time,  fill=dens))+
#   scale_fill_gradient2(low = "white" ,
#                        high = "red") +
#   main_theme